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Abstract. 

We formulate a theoretical model of the shear failure of a thin film tethered to a 
rigid substrate. The interface between film and substrate is modeled as a cohesive layer 
with randomly fluctuating shear strength/fracture energy. We demonstrate that, on 
scales large compared with the film thickness, the internal shear stresses acting on the 
interface can be approximated by a second-order gradient of the shear displacement 
across the interface. The model is used to study one-dimensional shear cracks, for 
which we evaluate the stress-dependent probability of nucleation of a critical crack. 
This is used to determine the interfacial shear strength as a function of film geometry 
and statistical properties of the interface. 



1. Introduction 



In this study we investigate the interfacial failure of thin films subject to shear loads. 
We consider situations where the interface with the substrate is disordered, such that 
the interfacial shear strength and fracture energy are random functions of position. 
While most of the mechanics literature has focused on the properties of interface 
cracks in systems where the properties of the interface are spatially homogeneous, 
failure of disordered interfaces with randomly fluctuating strength has been investigated 
extensively in the statistical physics community (see e.g. [U l2l El H] ) . Besides obvious 
applications to materials problems such as shear failure of coatings and shear-induced 
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delamination of thin films, the problem at hand has some interesting applications in 
geosciences [5], [6|, [7], where models similar to the one studied in this paper have been 
applied to the initiation of snow avalanches and landslides. 

In the present study we investigate theoretically how random variations of interface 
toughness affect the nucleation of interface cracks. We assume a one- dimensional 
geometry and evaluate the crack nucleation probability as a function of stress, geometry, 
and the statistical parameters characterizing the interfacial strength distribution. The 
theoretical arguments are validated by comparison with numerical simulations. 

2. Formulation of the model 

We consider a thin elastic film tethered to a rigid substrate. The interface with the 
substrate is modeled as a cohesive layer in the plane z = 0. The response of the 
interface to shear loads is characterized by a scalar stress-displacement relationship 
= 0) = r(-u) where axy is the shear stress at the interface, r is the interfacial 
shear strength, and u{x) = Wx{x,z = 0) is the shear displacement across the interface. 
The film is loaded in plane shear by spatially homogeneous tractions applied to its free 
surface, giving rise to a space- independent 'external' shear stresss 0"^^""". The maximum 
stress that can be supported by the interface is denoted as tm, and the failure energy 
is given by the integral 



where Uq denotes the characteristic displacement-to-failure. Stress-displacement 
relationships are schematically shown in Figure [H for the semi-brittle (full line) and 
brittle cases (dashed line). Structural disorder is introduced into the model in terms 
of statistical variations of tm which is considered a stochastic process with prescribed 
statistical properties that will be specified below. 

To analyze crack nucleation and propagation, we have to evaluate the internal 
stresses associated with the interfacial displacement field u{x). The elastic energy 
functional associated with a generic displacement vector field w{r) is given by 



where a = (2 — 2v) /(I — 2z/), and v is Poisson's number. Energy minimization gives the 
equilibrium equation 



We solve Eqn. [3] in Fourier space, imposing along the z = plane the boundary 
conditions Wx{x,z = 0) = m(x), Wy = = Q and making the crucial assumption that 
the film thickness D is much smaller than the characteristic length of variations in the 
displacement field: \kD\ << 1 where k is the wave vector modulus of any characteristic 
variations in w. The elastic energy can then be written as 




(1) 




(2) 




(3) 




(4) 
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Figure 1. Shear strength versus interfacial displacement (schematically); full line: 
semi-brittle behaviour, dashed line: brittle behaviour; dotted line: piecewise linear 
approximation to the semi-brittle case as used in our simulations. 



The total energy of the system is obtained by adding to the elastic energy H the 
work done by the shear stresses at the interface: 

fu{x) 



G{u) = -Jdx £^ (af^^^ - t{u)) du 



(5) 



Minimizing the total energy functional E{u) = G{u)+H{u) leads to the equilibrium 
condition 

where the gradient coefficient is given by / = 2D^{1 + a). Equations of this type have 
been studied by Aifantis and co-workers in the context of shear and slip bands in metal 
plasticity (see e.g. [S]). It may be noted that in these cases the constitutive relation 
contains a strain variable (shear strain or equivalent strain) in place of the displacement 
variable u. The mathematical structure is, however, the same as in the present problem. 



3. Energy of a shear crack and critical crack size 

We first consider interfaces with space independent fracture toughness. This serves to 
derive some relations and introduce notations that will prove useful in our subsequent 
analysis of the disordered system. We assume that a mode-II crack of finite length 
exists along the interface. Such a crack is characterized by a displacement field which, 
for X ^ oo , starts from a value Uq on the left stable branch of the t{u) curve, goes 
through a maximum Ui which without loss of generality we assume at x = 0, and then 
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reverts to Uq for x — oo. The derivation of the corresponding solutions of Eqn. 
has been discussed elsewhere |7] . In the limit of long cracks (small stresses) we find that 
the displacement profile is approximately parabolic, 

u{x) = II (7) 

where the parameter / defines the crack length. This equation describes the displacement 
profile as long as m 3> mq such that the shear strength t{u) ~ along the length of 
the crack. The parabolic profile is complemented by two boundary layers where the 
shear strength r(n) goes through the curve depicted in Figure [H For a long crack the 
contribution of these boundary layers can be neglected, and the energy is approximately 
given by 

^(0--^«^)' + 2^.,Z (8) 



This energy has a saddle point at the critical crack length Ic = y/2wfl/a. 



EXT 
xy ■ 

Cracks above this length are unstable under the load cr^^^ and lead to interface 
failure. The energy required to create a critical crack is Ec = iwflc/S, i.e., it decreases 
in inverse proportion with the applied stress. However, the barrier is always finite 
and, as emphasized by Arndt and Nattermann [4j for the case of bulk cracks, at 
physically meaningful stress levels the typical crack nucleation energies are many orders 
of magnitude above typical thermal energies. Hence, crack nucleation cannot occur as a 
result of thermally activated processes. However, energy barriers may be substantially 
reduced or even eliminated at certain locations due to the influence of structural disorder. 
In the following we analyse the conditions under which crack nucleation may occur 
spontaneously in a disordered medium. 



4. Shear crack nucleation at a disordered interface 

A disordered interface can be described by a randomly varying peak shear strength 
tm{x) or, equivalently, failure energy W{{x). In this case, the energy expression, Eqn. 
(IE]), modifies to 

m ^ {a^^^y + f_Wf{x)dx = [wf{x) - F{x)] dx (9) 

where F{x) = {xa^^'^^ / (2/) is the effective driving force acting on the edge of a crack 
of length X. 

A sufficient criterion for spontaneous nucleation of a propagating shear crack 
starting from the position x = is given by (dE/dl) < V/ or, equivalently, by 
F{1) > Wf{l) Wl. To estimate the probability for this to happen, we consider the case of 
short-range correlated disorder where the failure energy variations can be characterized 
by a finite correlation length ^. We split the interface into segments of length ^ and 
treat the shear strengths in each segment as independent, identically distributed random 
variables with probability distribution P{w{). The condition that the crack can advance 
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across the nth segment is F{nC,) > Wf{n^), the probabihty for this is P{F(n^)), and 
the crack nucleation probabihty is IVnuci = nn-P(-^(^0)- Taking the logarithm and 
reverting to continuous variables, we obtain 



lnW^nuci = ^/lnP(F(x))dx 



(10) 



To evaluate this probability we have to specify the probability distribution 
characterizing the variations of interfacial strength. We assume a WeibuU distribution 
with characteristic failure energy wq and modulus /3: 



P{w{) = 1 — exp 



Wo/ 



11 



The mean failure strength is then given by (wf) = woT{l + 1/(3) where T{x) denotes 
the Gamma function, and the coefficient of variation (GOV = variance/mean) of the 
distribution is calculated as GOV = [1(1 + 2//3)/{T{l + l/(3)f - 1]^/^. 

Using the driving force given above, the crack nucleation probability is then 
evaluated as 

1 



nucl 



lnP(F(a;))da; 



r(l + l//3) 



I ln[l - exp(-s2^)]ds 



2/(wf) 



xy 



9iP) 
(12) 



where s = a;//c is the ratio between the crack size and the size of a deterministic critical 
crack. Hence we obtain that 



nucl 



exp 



-9iP) 



21{wi) 

^'-'xy 



exp 



Er 



B-f eff 



(13) 



This has the structure of an Arrhenius equation where the activation energy is the 
energy of a critical crack in the system without disorder, and the place of temperature 
is taken by an effective 'disorder temperature' ksTes = 4:C,{wi) /3g{P). For large j3, we 
may approximate the integrand in Eqn. f|T2|) by Ins^^ for s < 1, and for s > 1. In 
physical terms, this means that only crack sizes up to the critical one significantly affect 
the nucleation probability, which is thus tantamount to the probability of forming a 
critical 'deterministic' crack. In this limit, we simply obtain that g{l3) = 2(3. With 
this approximation, the effective 'disorder temperature' is directly proportional to 
the variance of the WeibuU distribution times the correlation length of the strength 
variations. 

Nucleation of a propagating crack implies system failure. However, Eqn. (|T3ll 
cannot be directly interpreted as a system failure probability: While Eqn. f|T3l) is 
evaluated under the assumption that the crack starts from the origin, nucleation from 
any other position would produce a similar result. There are N = L/C, potential 
nucleation sites in a system of length L, and system failure occurs when a crack nucleates 
at the weakest of these. To evaluate the failure stress distribution, we use an argument 
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from extreme order statistics. Interpreted as a function of stress, W^iuciicr^^^) gives the 



probability that the system has, at any stress larger than cr^^^, produced a propagating 



crack starting from a given nucleation site. 1 — VFnuci(c"aiif ) is the probability that this 



has not yet happened, and [1 — W^^cii'^xy )] probability that the system as a 



whole is still intact at the stress cr^^'^- Hence we write the system failure probability as 



or 



Pia^^^) = 1 - exp [iVln(l - W^^,,)] (15) 



N 



xy 



For large systems, the nucleation probability at each site will be quite small, and we can 
therefore approximate ln(l — VFnuci) ~ W^nuci to obtain (up to logarithmic corrections in 
the exponent) 



P«^T) = l-exp , 



2I(wi) 



(16) 



This expression can be interpreted as the distribution of system failure stresses. 
Reciprocal failure stresses are Gumbel distributed, and the characteristic failure stress 
(e.g. the median of the failure stress distribution) scales like 



This stress decreases logarithmically with system size. It is interesting to compare 
this relationship with the deterministic failure stress of a system containing a crack of 
length I: cr^^^ = y/2wfT/l. It is evident that, in Eqn. (ITTl) . the correlation length of the 
disorder plays a very similar role to the crack size in case of a interface without disorder. 

5. Comparison with simulation results 

To test the above theoretical calculations, we use a lattice automaton technique where 
we evaluate the displacements at discrete sites Xi on a one-dimensional lattice with 
lattice constant C,- Accordingly, we replace Uxx in Eqn. ([6]) by the corresponding 
discrete second-order gradient. Furthermore we approximate the strength-displacement 
characteristics by a piecewise linear curve as shown by the dotted line in Figure [TJ 
Nondimensional variables are defined through 

EXT 

T = ^ S = f^ = U = ^,X = ^ (18) 

such that the average peak strength and fracture energy are by definition equal to 1. In 
these coordinates the equilibrium equation reads 

J^ + S-T<0 (19) 
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where the interaction constant J is given by J = {Iuq)/ (tm^^). The length of a critical 



and the failure strength distribution of the disordered system, Eq. ( fT6l) . reads 



In our simulations, random values for the local peak strengths are generated as 
Weibull distributed random fields with average 1, spatial correlation length 1 and 
WeibuU parameter /3. Note that the nondimensional value of the parameter /?, together 
with the interaction constant J and the system size N , are the only independent 
parameters of the problem. 

The system is slowly loaded by increasing the external stress S from zero in small 
steps until sites become unstable as the local (external plus internal) stress exceeds 
the local shear strength. The displacement at all unstable sites is then increased by a 
small amount Af/. New internal stresses are computed for all sites and it is checked 
again where the sum of external and internal stresses exceeds the local strength. The 
displacement at the now unstable sites is further increased, and the process is repeated 
until the system has reached a new stable configuration. Then the external stress is 
increased again and so on until the system has failed completely {U > 2 for all sites). 
The corresponding critical stress is recorded, and the procedure is repeated for different 
realizations. 

Figure [2] shows a failure stress distribution determined for J = /? = 1 in a system 
of size = 100. Fitting the functional form of Eq. ( l20i) by using g{(3)\/2J as a 
fit parameter yields g{(3) = 2.12 as compared to g{P) = 1.64 obtained by direct 
computation. The discrepancy is acceptable in view of the approximations made in 
the derivation of the failure stress distribution. Using the same parameter, the size 
dependence of the mean failure stress is well reproduced (inset in Figure El). 

Figure [3] shows the dependence of the mean failure strength on the value of the 
interaction coefficient and on the degree of disorder. As seen in the main graph, the 
decrease of strength with increasing disorder (increasing CoV/decreasing (3) is in the 
large-disorder regime well described by the theory In the regime of small disorder, the 
theory predicts a diverging strength whereas the simulation for zero disorder yields a 
finite strength of 1. This is due to the finite value of the maximum stress that can be 
sustained by the interface; at this point the system fails by homogeneous yielding rather 
than by crack nucleation [for ideally brittle behavior the maximum stress would be 
infinite]. The effect can possibly be taken into account by accounting, in the derivation 
of the crack properties, for a finite process zone size as discussed by Bazant [9]. 

According to the theory, the characteristic failure stress is expected to scale as 
S oc g{l3)\/2J / ln(N). The inset shows that failure stresses determined for different 
system sizes and Weibull exponents indeed collapse when scaled by g{P)/\nN and 
plotted vs. J, even though the fit function yields an exponent of 0.38 rather than 



crack in a homogeneous system is in non-dimensional coordinates given by 



V2J/S, 



P(S) = 1 - exp [exp (-(7(/3)Le(S) + In(Ar))] 




(20) 
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Figure 2. Main grapli: failure stress distribution, parameters: J — (3 — \, N ~ 100, 
1000 simulations; full line: theoretical fit using Eq. ^] the fit yields g(/3)%/2J = 3.00. 
Inset: Size dependence of the mean failure stress, parameters J — [3 ~ \] full line: 
E = 3/(ln(iV) + 3). 




Figure 3. Main graph: mean failure strength vs. coefficient of variation, parameters: 
J = l,iV = 1000; full line: theoretical fit (plot of 0.16g(/3) vs. GOV). Inset: mean 
failure strength, normalized by g{[3)/\n{N), vs interaction coefficient J; full line: fit 
curve 1.15J°-^^. 
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0.5. This discrepancy is attributed to logarithmic correction terms neglected in our 
derivation of the failure strength distribution. 

6. Conclusions 

We have derived an analytical expression for the nucleation probability of one- 
dimensional cracks in a disordered interface between a thin film and a rigid substrate. 
The method we used is a variation of the method used by Arndt and Nattermann |1] for 
deriving the nucleation probability of bulk cracks. Accordingly the obtained relations 
have strong structural similarities with those found in P]. A main difference between 
the present problem and the one treated by Arndt and Nattermann resides, however, 
in the nature of the elastic interactions. Elasticity is essentially non-local for a bulk 
crack whereas, in the present case, the small thickness of the film leads to an effectively 
local stress re-distribution. From this point of view, our model may be considered a 
continuum version of local-load-sharing (LLS) fibre bundle models. Failure probability 
distributions for such models have been derived by several authors (see e.g [IDj and 
references therein). For certain types of disorder, such as a binary strength distribution, 
exact results can be obtained which confirm the existence of a logarithmic size effect. The 
present treatment, which uses extreme order statistics to obtain an expression for the 
size-dependent probability of system failure, is much more approximative. However, it 
has the advantage that the methods used can be applied to a broad class of local strength 
distributions with support in the interval [0, oo]. Applying the present method to other 
local strength distributions (e.g. log-normal) leads to a similar structure of the failure 
probability distribution. The specific functional form of the distribution influences only 
the calculation of the function g which characterizes the relation between the failure 
strength and the parameter governing the distribution width. 

The present investigation focuses on one-dimensional cracks. This is a major 
restriction since in one dimension the crack shape does not depend on the disorder. 
In two dimensions, on the other hand, the crack shape may adjust to the strength 
fluctuations. In this case, the configuration space for crack nucleation becomes 
much more complicated as the problem can no longer be treated in terms of the 
crack length only, which in the present work serves as a single reaction coordinate. 
Instead, an investigation of crack nucleation and propagation in two dimensions requires 
consideration of the interplay between disorder and crack front elasticity, and a 
treatment of the ensuing crack front roughening. These aspects will be considered 
in a later publication. 
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